% 
% posBS = [0, 0]; % 基站在二维坐标系中的位置
% 
% posUE = [10; 0]; % 近场用户在二维坐标系中的位置
% 
% posEV = [0; 10]; % 窃听者在二维坐标系中的位置 
% 
% D = 0.3; % 整个均匀平面阵列的孔径
% 
% lambda = 1e-3;  % 毫米波波长30Ghz
% Nd = floor(2 * D / lambda); % 均匀平面阵列水平方向阵元数
% Ne = Nd; % 均匀平面阵列垂直方向阵元数
% spacing = D / Nd;  % 阵列阵元之间的间距
% 
% 
% % 生成均匀平面阵列的阵元坐标
% posBSelements = zeros(Nd, Ne, 2);
% for i = 1:Nd
%     for j = 1:Ne
%         x = (i - 1) * spacing - D / 2; % X坐标
%         y = (j - 1) * spacing - D / 2; % Y坐标
%         posBSelements(i, j, 1) = x;
%         posBSelements(i, j, 2) = y;
%     end
% end
% 
% 
% 
% 
% posBSelements = zeros(Nd, Ne, 2);
% size(posBSelements(1,1,:))
% 
% 
% alpha_c = 1.2;
% beta_c = 827.67;
% Pmax = 1;


%% 信道系数的优化

% 定义参数和变量
D = 0.3; % 整个均匀平面阵列的孔径
lambda = 1e-2; % 波长0.01m，单位：米   30Ghz
N_d = floor(2 * D / lambda); % 垂直方向的阵元数  600
N_e = floor(2 * D / lambda); % 水平方向的阵元数  600
spacing = D / N_d;
M = 2; % 用户数
b = 2; % Boresight gain参数
x_b = 5; % 用户 Bob 的 x 坐标
y_b = 5; % 用户 Bob 的 y 坐标
z_b = 0; % 用户 Bob 的 z 坐标
x_w = 3; % 窃听者 Willie 的 x 坐标
y_w = 8; % 窃听者 Willie 的 y 坐标
z_w = 0; % 窃听者 Willie 的 z 坐标
k = 2 * pi / lambda; % 波数
alpha_c = 1.2;
beta_c = 827.67;


% 生成信道系数
% rho = rand(N_d, N_e); % 元件位置
A_b = zeros(N_d, N_e);
A_w = zeros(N_d, N_e);
s = randn(N_d, N_e); % 模拟发送信号

for i = 1:N_d
    for l = 1:N_e
        % rho(i, l) = 
        p_i_l = [(l - 1) * spacing - D / 2, (i - 1) * spacing - D / 2, 0];
        % 计算用户 Bob 的信道系数
        p_b = [x_b, y_b, z_b];
        Theta_i_l_b = calculateTheta(p_i_l, p_b);
        distance_b_il = norm(p_b - p_i_l); % 计算距离 |p_m - p_i_l|
        phase_b_il = exp(-1i * k * distance_b_il); % 计算相位 e^{-j*k*|p_m - p_i_l|}
        A_b(i, l) = sqrt(F(Theta_i_l_b)) * (lambda / (4 * pi * norm(p_b - p_i_l))) * phase_b_il;

        % 计算窃听者 Willie 的信道系数
        p_w = [x_w, y_w, z_w];
        Theta_i_l_w = calculateTheta(p_i_l, p_w);
        distance_w_il = norm(p_w - p_i_l); % 计算距离 |p_m - p_i_l|
        phase_w_il = exp(-1i * k * distance_w_il); % 计算相位 e^{-j*k*|p_m - p_i_l|}
        A_w(i, l) = sqrt(F(Theta_i_l_w)) * (lambda / (4 * pi * norm(p_w - p_i_l))) * phase_w_il;
    end
end

% 计算接收信号
a_b = A_b(:); % 将 A_b 转化为列向量
a_w = A_w(:); % 将 A_w 转化为列向量


%% DMA的发送信号
N = N_d * N_e; % 总元件数量
sigma_b2 = 1; % 噪声方差

% 生成随机的数字前向编码向量
w_b = (randn(N_d, 1) + 1i * randn(N_d, 1));
w_w = (randn(N_d, 1) + 1i * randn(N_d, 1));


% 初始化 H 矩阵和 Q 矩阵
H = zeros(N, N);
Q = zeros(N, N_d);

% 初始化示例的元件相位和位置


phi = linspace(0, 2 * pi, N_d); % 元件响应角   待优化的变量

% 填充 H 矩阵和 Q 矩阵
for i = 1:N_d
    for l = 1:N_e
        % 计算元件相位
        
        h_i_l = exp(  ((l-1)*spacing-D/2)  *  (alpha_c + 1i*beta_c)  );
        
        % 计算可调响应 q_i_l
        q_i_l = (j + exp(1i * phi(i))) / 2;
        
        % 将 h_i_l 填充到 H 矩阵对角线
        H((i - 1) * N_e + l, (i - 1) * N_e + l) = h_i_l;
        
        % 将 q_i_l 填充到 Q 矩阵
        Q((i - 1) * N_e + l, i) = q_i_l;
    end
end

% 生成示例输入信号 x_b
x_b = ( randn(1, 1) + 1i * randn(1, 1) ) / sqrt(2);
x_w = ( randn(1, 1) + 1i * randn(1, 1) ) / sqrt(2);

s = H * Q * (w_b * x_b + w_w * x_w);

% % 计算可达速率
% a_m = ones(1, N); % 示例 a_m 向量
% w_j = ones(N, 1); % 示例 w_j 向量
% R_m = log2(1 + abs(a_m * (H * Q * w_j))^2 / (sum(abs(a_m * (H * Q * w_j))^2) + sigma2));
% 
% disp(['可达速率 R_m = ' num2str(R_m)]);

s = s(:); % 将发送信号转化为列向量
n_b = sqrt(0.5) * (randn(1, 1) + 1i * randn(1, 1)); % 用户 Bob 的噪声
n_w = sqrt(0.5) * (randn(1, 1) + 1i * randn(1, 1)); % 窃听者 Willie 的噪声
r_b = a_b' * s + n_b; % 用户 Bob 的接收信号  
r_w = a_w' * s + n_w; % 窃听者 Willie 的接收信号

disp('用户 Bob 的接收信号:');
disp(r_b);

disp('窃听者 Willie 的接收信号:');
disp(r_w);


%% 待优化变量： Pb Q w_b w_w
Pbmax = 10;
Pb = 1; % 只剩Pb没初始化了， variable Bob的发射功率待优化
epsilon = 0.01;

% 生成一个在[0, Pwmax]范围内的均匀分布的Pw
Pwmax = 2;
Pw = Pwmax * rand;


% Bob处接收到的信噪比
gamma_b = Pb * norm(a_b' * H * Q * w_b)^2 / ( Pw * norm(a_b' * H * Q * w_w)^2 +  sigma_b2) ;  % aim to max this
disp('Bob 的信噪比:');
disp(gamma_b);

gamma_w = Pwmax/2 * norm(a_w' * H * Q * w_w)^2 / ( Pw * norm(a_w' * H * Q * w_w)^2 +  sigma_b2) ;  % aim to max this
disp('Bob 的信噪比:');
disp(gamma_b);

% 
% 
% % constrains
% % Decision Error Probility of Willie
% Pr_b = 1 - Pb * norm(a_w' * H * Q * w_b)^2 / ( Pwmax * norm(a_w' * H * Q * w_w)^2 );
% 
% Pr_b >=  epsilon;
% 
% Pb <= Pbmax;
% 
% norm(w_b) == 1;
% 
% norm(w_w) == 1;
% 
% % Q的约束
% 









